On-site number statistics of ultracold lattice bosons 



o 
o 

(N 
D 

m 
(N 



Barbara Capogrosso-Sansone/ Evgeny Kozik,^ Nikolay Prokof'ev,^'^ and Boris Svistunov^'^ 

''Department of Physics, University of Massachusetts, Amherst, MA 01003 
^Russian Research Center "Kurchatov Institute", 123182 Moscow, Russia 

We study on-site occupation number fluctuations in a system of interacting bosons in an optical 
lattice. The ground-state distribution is obtained analytically in the limiting cases of strong and 
weak interaction, and by means of exact Monte Carlo simulations in the strongly correlated regime. 
As the interaction is increased, the distribution evolves from Poissonian in the non-interacting gas 
to a sharply peaked distribution in the Mott-insulator (MI) regime. In the special case of large 
occupation numbers, we demonstrate analytically and check numerically that there exists a wide 
interval of interaction strength, in which the on-site number fluctuations remain Gaussian and are 
gradually squeezed until they are of order unity near the superfluid (SF)-MI transition. Recently, 
the on-site number statistics were studied experimentally in a wide range of lattice potential depths 
[Phys. Rev. Lett. 96, 090401 (2006)]. In our simulations, we are able to directly reproduce 
experimental conditions using temperature as the only free parameter. Pronounced temperature 
dependence suggests that measurements of on-site atom number fluctuations can be employed as a 
reliable method of thermometry in both SF and MI regimes. 
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I. INTRODUCTION 

Experiments on ultracold atoms trapped by an opti- 
cal potential Q, Q, IE Hi Uli offer a unique possibility 
to explore fundamental properties of strongly correlated 
quantum many-body systems allowing virtually unlim- 
ited control over the microscopic Hamiltonian parame- 
ters (see, e.g., Q and references therein). System flexi- 
bility along with relatively long decoherence times puts 
it among top candidates for implementation of quantum 
information algorithms 6]. In atomic interferometry 0, 
ultracold gases in the strongly correlated regime allow 
to achieve accuracies below the standard shot noise limit 
[E 1^ ' A new exciting application, "atomtronics" , is sug- 
gested by a remarkable analogy between the physics of 
ultracold atoms in optical lattices and that of electrons in 
crystals 'lO^. With the current experimental technique, it 
seems plausible to produce such basic atomtronic devises 
as diodes and bipolar junction transistors, which serve as 
building blocks for amplifiers and logic gates. 

On the fundamental physics side, strongly interacting 
lattice bosons provide insight into the nature of quantum 
phase transitions. In particular, these systems are an 
accurate experimental realization of the Bose-Hubbard 
model [13, 



H 



<hj> 



a] a. 



1) 



E 



(1) 



where hi = a^a^ is the number operator on site i; a| and 
Ui respectively create and annihilate bosons on lattice 
sites, and < i,j > denotes the sum over nearest neigh- 
bors. The first term describes tunneling between neigh- 
boring potential wells of the optical lattice, the second 
term is the effective repulsion within a well, while the 
last term is due to an additional smooth space-varying 
potential, such as, e.g., a magnetic trap. This system 



exhibits a transition between superfluid (SF) and Mott- 
insulator (MI) groundstates governed by the competition 
of atom mobility and interatomic interaction p^ . The 
SF-MI phase transition is a topic of intense current re- 
search both theoretically 0, Q, 0| and experimentally 
[TgL im IT^ . In typical experiments, a prepared Bose- 
Einstein condensate of ultracold atoms is driven to the 
strongly correlated regime by means of its adiabatic load- 
ing into a periodic optical potential (optical lattice) in- 
duced by the a.c. Stark effect of interfering laser beams. 
The mobility of atoms, namely the hopping t, and their 
interaction U are controlled by the depth of the optical 
potential, i.e. by the laser intensity. In relatively shallow 
potentials, atoms are delocalized over the entire lattice 
giving rise to long-range coherence and thus to su- 
perfluid behavior. 

If the lattice filling is commensurate, i.e. there is on 
average an integer number of atoms per lattice site, in- 
creasing the lattice depth brings the system across the 
phase transition to the MI state, which is characterized 
by zero compressibility and a gap in the spectrum of el- 
ementary excitations. Here, the key observable is the 
atom interference pattern obtained upon releasing the 
atoms and letting the atom cloud expand for a transient 
time of flight [J . Phase correlations between lattice sites 
result in pronounced interference peaks smeared by the 
finite correlation length in the MI regime. At the same 
time phase coherence is fundamentally connected with 
the statistics of atom number fluctuations on lattice sites. 
In particular, one expects the on-site atom number to 
behave as a canonically conjugate variable with respect 
to the phase field and, therefore, experience suppressed 
fiuctuations in the MI regime, analogous to the num- 
ber squeezed states with sub-Poissonian number fluctu- 
ations 2Ql widely studied in quantum optics (see, e.g., 
[2lj|). in practice, number squeezed states are impor- 
tant for high- precision atomic interferometry Q, where 
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their use can potentially lead to sensitivities limited only 
by the Heisenberg uncertainty principle [sL I^j23|. and 
for atom-based quantum computing techniques Q , where 
unwanted number fluctuations necessitate correction pro- 
cedures in operation of quantum gates. 

Until recently, the number distribution was not mea- 
sured directly. The situation changed with the develop- 
ment of the spin-oscillation technique 3] , which is sensi- 
tive to the number of atom pairs and works at arbitrary 
lattice depths [ITj . and, most recently, the microwave 
spectroscopy using atomic clock shifts 4J. In Ref. [l^ . 
Gerbier et al. observed a drastic change of atom number 
statistics as the system of ^^Rb atoms was driven through 
the SF-MI transition. On the theoretical side, apart from 
the recent mean- field calculation |23l |. a comprehensive 
study of atom number fluctuations in the strongly corre- 
lated regime is still missing. 

In the present work, we attempt to close this gap by 
tackling the problem both analytically and numerically. 
In section^ we focus on the academic case of a homoge- 
neous square lattice in the thermodynamic limit and in 
the limit of zero temperature in one-(lD), two-(2D) and 
three- (3D) dimensions. At ?7 = 0, the ground state of an 
ideal Bose gas is a Bose-Einstein condensate with char- 
acteristic Poisson distribution of number fluctuations. In 
subsection III Al we consider the limit of weak interaction 
uU/t <^ 1, where v is the filling factor. This parame- 
ter region corresponds to perturbative squeezing of the 
Poisson distribution, which, as a function of vU /t, quali- 
tatively depends on the space dimension d. We consider 
the case of large filling factors ^ 1 separately in sub- 
section]^^ because it is qualitatively different from that 
oi V ~ 1. Indeed, the quantum nature of the SF-MI 
transition implies that the number fluctuations in the 
vicinity of the transition must be of order unity. At the 
same time, sX U — Q the variance of the on-site num- 
ber distribution is, — v ^ \. Therefore, there exists 
an extensive range of interactions (defined by the condi- 
tion U/t < z/), in which the system remains superfluid, 
but its on-site number distribution is drastically squeezed 
before the SF-MI can take place. We show that, at 
U/t <^ I', the on-site number statistics are Gaussian and 
derive the variance of the distribution, which scales 
as oc ^ vtjU at 1/:^ <C U jt ^ z/ in all dimensions. In 
subsection III CI we suggest an expression that interpo- 
lates P{ri) between the limiting cases of small interaction 
and large occupation numbers, which is found to properly 
describe P(n) up to C//t of the order of the critical value 
{U /t)c- The strong coupling limit, vt/U <C 1, at integer 
filling is considered in subsection III Dl In this limit, the 
system is in the MI regime and the on-site number dis- 
tribution is governed by rare particle-hole fluctuations. 

We study the distribution in the strongly correlated 
regime connecting the limiting cases by means of a di- 
rect numeric simulation of the model ^ &t v — 1 in 
subsection III El The distribution of the on-site occu- 
pation number is a local property, revealing no critical 
features at the transition. However, the strongly corre- 



lated region is responsible for the crossover that changes 
the statistics qualitatively. As the interaction strength is 
increased, we observe a gradual squeezing of the on-site 
number distribution and the emergence of the symmetry 
between particle- and hole-like fluctuations, characteris- 
tic of a MI. In this section, we also present numerical 
data for the case of large filling factors and demonstrate 
that the analysis of subsection III Bl is applicable already 
a,t V = b. 

The worm algorithm quantum Monte Carlo (MC) tech- 
nique (24] allows us to simulate system sizes that are 
currently realized in experiments without any approxi- 
mations, including the particle number. The results of 
a direct numeric simulation of the experimental setup 
of Ref. |l/3| are presented in section IIIII With the lat- 
tice parameters fixed by the experiment, we are left with 
temperature as the only free parameter. Due to a smooth 
confining potential present on top of the optical lattice, 
the number distribution is not homogeneous. We focus 
on an integral characteristic of the number distribution, 
namely, the fraction of atoms found on lattice sites with 
occupation n, which can be systematically measured ex- 
perimentally. This quantity has a pronounced tempera- 
ture dependence in both SF and MI regimes. 

The problem of thermometry in optical lattices, espe- 
cially in the MI regime, is a long standing one. The abil- 
ity to control the temperature is of crucial importance 
for applications that rely on the peculiar properties of 
a MI state. At T = fluctuations are of quantum na- 
ture and can be efficiently controlled externally through 
the lattice parameters, whereas temperatures compara- 
ble to the Mott excitation gap destroy the insulating 
state by activating particle-hole excitations. At the mo- 
ment, there are no experimental techniques to measure 
the temperature of a strongly interacting system. Unlike 
in the weakly interacting regime, where the temperature 
can be straightforwardly extracted from the momentum 
distribution (e.g., from the interference pattern of mat- 
ter waves or from the condensate fraction observed after 
the trap is released and the gas expands freely), in the 
strongly correlated regime, both temperature and inter- 
atomic interaction are responsible for filling the higher 
momentum states making standard absorption imaging 
techniques inapplicable. 

The idea of using occupation number distributions to 
estimate the temperature was explored in Ref. |25j , where 
in was argued that the temperature dependence of the 
total number of pairs and their spatial distribution (in 
traps) provides a sensitive method of thermometry deeply 
in the MI phase at energies smaller than the interatomic 
interaction, but much larger than the effective hopping 
between the sites. [Recently, it has become possible to 
directly sample spatially-resolved number distributions 
by spin changing collisions |l6j |. microwave spectroscopy 
4], and the scanning electron microscope '26'| promis- 
ing a complete practical realization of this method.] In 
this paper, we perform thermometry of the system in all 
strongly correlated regimes by comparing experimental 
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data and numerical results for the statistics of occupa- 
tion numbers. More specifically, we compare numerically 
calculated fraction of pairs {n — 2) with that measured 
in Ref. 0| across the SF-MI transition estimating the 
range of experimental temperatures. The accuracy of 
this method is mainly limited by the error bars of the 
experimental data and by the range of applicability of 
the Bose-Hubbard model. 



II. 



HOMOGENEOUS LATTICE 



In this section, we assume that there is no space- 
varying potential on top of the optical lattice, and set 
Ei = eQ. Let Ns be the number of lattice sites and N be 
the total number of particles. The goal of this section is 
to obtain the ground state probability P„ to detect n par- 
ticles on a given lattice site in the limit oi Ns, N oo, 
at a fixed filling factor v = N/Ng. Mathematically, P„ 
can be defined as 



Pn{U/t) 



E I 



{m = ■n,{n,^i}\ \'^u/t)\ 



(2) 



where |5'c//t) is the many-body ground state wavefunc- 
tion, and |{?T.i}) are Fock states. 



A. Weak coupling limit 

At sufficiently small U/t, the relevant representation 
of Eq. is obtained by the diagonalization of the ki- 
netic energy term with the canonical transformation a.; ~ 



N. 



-1/2 



X^k exp(i27r kr^/i) (periodic boundary condi- 



tions are assumed), where L = Ng^'^ is the linear system 
size, Ti and k are respectively the position of the site 
i and a quasi-momentum in d dimensions with integer- 
valued components, — (L — l)/2 < r^^, fc^ < {L — l)/2, 
11= 1, ... ,0?. The result is 
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ki+k2=k3+k4 



with ek = 2tX;^=i [1 " cos(27rfc^/L)]. At U = the 
ground state of the Hamiltonian Q is a pure Bose- 
Einstein condensate, which can be expressed as a co- 
herent state, |*((7/t)=o) = exp(VlV6j - iV/2)|0). Then, 
transforming the wavefunction back to the on-site repre- 
sentation in terms of {oi} yields the Poisson distribution 
for the probability to find n particles on a given site: 



n 



(4) 



At a finite, but small, U we employ the standard Bo- 
goliubov method of separating the system into the 
classical-field condensate part and non-condensate parti- 
cles interacting with it, omitting the terms of the third 



and forth order with respect to the non-condensate op- 
erators. In this approximation, the Hamiltonian (O is 
reduced to a bilinear in fej. and form and diagonal- 

ized by the canonical transformation Ck = Uk^k + I'k^lk' 
where 

wk= [(ek + t^C/)/2L^k + l/2]'/', 
vk = [(£k + iyU)/2ujk - l/2f^ , 

u;^=[el + 2iyeuUY^\ (5) 

The ground-state wavefunction is then obtained from the 
equation Ck l^'cz/t) = for all k 7^ and has the form 



\^u/t) = Cexp 



2 ^ — ' Mk 

kT^O 



|0), (6) 



where C is the normalization factor and Nq = N — 
X^kT^o^^k^k) is the number of condensate particles. 

Now we can express l^^jj/t) in terms of the on-site op- 
erators. 



I* 



u/t) 



Cexp 



2 



with i/Q = Nq/Ns and 



27V, ^ Mk 



|0), 
(7) 

(8) 



In this form, the wavefunction can be used to obtain 
the on-site number distribution in the whole range of 
U/t <C by a straightforward application of Eq. ((SJ. 

We derive a closed-form expression for P(n) in the lim- 
iting case of 



vU/t < 1, 



(9) 



which corresponds to the range of sufhciently weak 
squeezing of P(n) allowing us to consider only the first 
correction to the Poisson distribution in the leading 
power of a. Mathematically, Eq. © implies that Sij <C 
1. If, in addition, Sij is short-range, i.e. it decays at 



distances Ir,; 



1, which implies that the leading 



correction is insensitive not only to the system size, but 
also to the value of the healing length cx then 
we can expand the exponential in Eq. in powers of 
Sij. Rather straightforward but lengthy algebra yields 
the distribution in the form 



p - p(0) 



f \(a) 



[^(0) _ 2P("), + p("),] 



(10) 



where A is an interaction-dependent squeezing parameter 
and Pn^ is given by Eq. assuming Pn^ = for n < 0. 
Identically, 



A(«) 



n ^ {n — 



(11) 
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Let us postpone writing an explicit expression for 
A and discuss the underlying assumptions leading to 
Eqs. H10 |I . H11|I . It turns out that Sij is local only in 3D, 
where the main contribution to the integral in Eq. (jS)) 
comes from large momenta close to the edge of the Bril- 
louin zone, |k| ~ L/2. In 2D, the integral has a logarith- 
mic divergency at low momenta in the limit a —f 0, while, 
in ID, the dominant contribution to Sij comes from small 
momenta, meaning that in both ID and 2D cases we can 
not rely on the perturbative expansion of the wave func- 
tion. Nevertheless, Eqs. ((Tmi . lfTTl) are still valid in ID 
and 2D, since the functional form of P{n) must be the 
same in all dimensions. To prove this statement, we note 
that P{n) is unambiguously determined by its character- 
istic function x{t) — (exp(j i a|ai)), which can be used 
to generate all moments of P{n). The function x{t) can 
be calculated explicitly as a series expansion in powers 
of {italui). Since averages of bosonic quasiparticle op- 
erators in the Bogoliubov theory obey Wick's theorem, 
each term in the series is a function of ^ = (ajai) and 



(aitti). Therefore, x(t) = x(^;C;C'); meaning that 



all physical parameters, including the space dimension, 
enter P{n) only through ( and In Eqs. ifTUI) . lfTT|l . this 
combination determines A(a). 

To determine A (a) in all dimensions were note that it 
is directly related to the dispersion of P{n) in Eq. 



(12) 



On the other hand, = {nf) — = (ala^ala^) — 
can be calculated in a standard way by replacing 
with their expressions in terms of the Bogoliubov modes 
and the classical-field condensate part, i.e. = JTi^ -|- 



Ws 



'1/2 



Sk^^oI^fcCk-WfccL^] exp(i27rkri/L). Then, an ap- 
plication of the Wick's theorem along with (ckcj^,) = (5^ k' 
and (ckCk') = gives 



^ tk 



£k 



(13) 



Strictly speaking, Eq. H13|) is valid as long as the num- 
ber of non-condensed particles is small (the ID case is 
special in this regard and is further discussed below), 
{N~No)/N < 1, which implies U/t < vi^-d^y^ ^^^-^^ 
V and U/t <^ V &% large v. In the latter case, the dis- 
tribution change can be non-perturbative since A ~ 1 is 
typical in this parameter regime, which will be discussed 
in more detail in the next subsection. 

Thus, we arrive at the following result 



Ns rv„ ^ ^k^ 



(14) 



k/O 



The asymptotic behavior of A(a 0) qualitatively de- 
pends on the dimensionality. In ID, the main contribu- 
tion to the integral H14|l comes from low momenta result- 
ing in 



A 



(a) 



V2 



la. 



(15) 



In 3D, there is no low- momentum singularity at a — > 0, 
and the asymptotic expression is linear in a: 



B 



(16) 



B 



dzi dz2 dz3 



oJoJo Eu=i(l - cosz^) 



15.673. 



(17) 



In 2D, the low- momentum singularity is logarithmic 



A 



(d=2) 



(a) 



ln(C/a) 



47r 



C « 23.54 



(18) 



A comment is in order here concerning the derivation 
procedure for the ID case. Formally, the condensate 
fraction is zero in the macroscopic limit even at T = 0, 
while the derivation is based on the assumption that al- 
most all the particles are Bose condensed. Nevertheless, 
our final results for the probabilities are valid even in 
ID, and the generic reasoning — based on the notion of 
quasicondensate — leading to this conclusion is as follows 
28]. The quasicondensate correlation properties charac- 
teristic of a weekly interacting ID gas at T = imply 
two different correlation radii, Tc and Rc, Rc ^ re- Here 
Tc defines the length scale upon which the system can 
be considered as macroscopic, while Rc is the length at 
which (quantum) fiuctuations of phase are of order unity. 
If the system size L is such that 



<^L<^ Rc 



(19) 



then the system is macroscopic with respect to all lo- 
cal properties, while still featuring a genuine condensate. 
(In a ID weakly interacting system at T = the density 
of this condensate is close to the total density of parti- 
cles.) Hence, for all local properties, including the ones 
discussed in the present paper, one can assume, without 
loss of generality, that the system size is finite and satis- 
fies the condition (|19|) . It should be emphasized that this 
assumption is implicit and is used exclusively to justify 
the derivation procedure. Otherwise, it does not lead to 
any explicit dependence of final answers on L. Indeed, 
the first inequality in Eq. (|19|l guarantees that all discrete 
sums over momenta can be replaced with integrals, and, 
since the integrals are convergent, the answer is indepen- 
dent of L. 



B. Large occupation number limit 

At i/ 1 there exists a wide {U/t <C i^) superfluid re- 
gion, in which large number fluctuations given by Eq. Q 
are gradually suppressed by the interaction until they be- 
come of order unity at the SF-MI transition. This physi- 
cally appealing regime is not captured by Eqs. H10() . H11|I . 
since they are applicable only for U/t <C l/j^. 

At v ^ 1 and U/vt <C 1, the number distribution is 
easily obtained due to the fact that the typical values 
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of the occupation number fluctuations, rji = ajai — v, 
are large 1 <C \rji\ <C v. Thus, the transformation 
ai — + r\i exp(i$i), where and $j are canonically 
conjugate Hcrmitian operators, along with rji/v <^ 1 re- 
duces the Bose-Hubbard Hamiltonian |^ to 



cos($i — $j 



+ ^T.^l (20) 



where < $i < 27r. In fact, in this form, Eq. ll^ is 
applicable at any iji <^ v including the description of the 
SF-MI transition and the MI phase. At sufficiently large 
interaction, namely U/t ^ l/i/, the oc rjf/i/'^ term in 
Eq. (|2()|l can be neglected, and the Hamiltonian coincides 
with the well-studied quantum rotor model (QRM). 

In this work, we are interested in the U/t <^ v limit 
of model Eq. 1)20(1 . i.e. when phase fluctuations between 
the nearest-neighbor sites are small |$i — $j | ^ 1 and 
the number fluctuations are large. This allows us to con- 
sider rji as a continuous variable and formally redefine 
the domain of 77^, $i as —00 < rii,^i < 00. The result is 
a harmonic approximation of the Hamiltonian (|20|) . 



H 



U 



-if- 



From its functional (quadratic) form we immediately 
conclude that the distribution of rji, which is Gaussian 
at f7 = (i^ ^ 00 limit of Eq. I^J), remains Gaussian 
in a wide range of coupling strength — all the way to the 
vicinity of the SF-MI transition where ct^ ~ 1 and the 
model (I21|l breaks down. The proof is as follows. Since 
the Hamiltonian 1(2 l|l is bilinear in {77^,$^}, all averages 
are subject to the Wick's theorem. Therefore, the charac- 
teristic function of the distribution W{rii), which is given 
by the integral /_ exp(ikr]i)W {r]i) dr/i, is a Gaussian 

exp(-fc2aV2). 

The only parameter of the distribution, cr^, is given by 



£k 



1 + 2(4ck) 



(22) 



where Ck and are the creation and annihilation oper- 
ators of the normal modes of the Hamiltonian l(21() with 
frequencies Wk given by Eq. |SJl. The Hamiltonian is di- 



agonalized by the transformation 



E 

k 



Wk 



2[/ + £k/j^ 



^E 

k 

Xki 



^k 
2i/ek 

1 



XkjCk-Xki4 



exp(i27r kr^/L) 



(23) 



Here, we restrict ourselves to ground state properties 
only and thus set (C]^Ck) = 0, which leads exactly to 
Eq. 1(1^ . This is hardly surprising, since the model 
is equivalent to the Bogoliubov approximation of the 
Hamiltonian in the limit of large v, and thus, we could 
formally demonstrate that P{n) is Gaussian a,t v ^ 1 al- 
ready in the framework of section FlI AI However, the hy- 
drodynamic approach chosen in this section seems more 
natural and physically transparent when dealing with a 
dense system. 

Let us examine properties of the distribution variance 
in more detail. An explicit expression for reads 



~d 




E^(l-cosa:^) 
vU/t + X!u(l ~ COSXp) 



dxi ■ ■ ■ dxd 



(24) 

In the limit of a = vU /t 0, this expression reduces to 
Eqs. ((15(l - l(18() . In the opposite limit, we find 



(25) 



where h = 2\/2, h ~ 13.373, and I3 « 52.348. The 
latter formula allows to estimate the range of U/t at 
(^-'^) which Eq. is apphcable. The condition > 1 

gives U/t <C V, consistent with the applicability range 
of Eq. 



C. Interpolation formula in the SF regime 

In sections III AI and IIIBI we have derived asymptoti- 
cally exact solutions describing week squeezing (A <C 1) 
and strong squeezing (0<A<l)ati/;»l respectively. 
The two limits overlap, which suggests that one can write 
a single interpolation formula to capture both limits cor- 
rectly. This formula is also expected to predict P{n) cor- 
rectly in a broader region of U/t even at v ^ 1. Formally, 
we have to find a function P{n, v, A) such that (i) it coin- 
cides with Ea. pi(l in the limit of A <C 1, (ii) it becomes 
Gaussian as a function of n at ^ 1 and < A < 1 
with the average v and variance = 1^(1 — A), and (iii) 
it is analytic with respect to u and A. The solution is 
not unique, and we simply suggest one which satisfies 
the above mentioned criteria 



Pin,iy, A) = 
exp 



cP^°Hn) 

\{n — v) + {n — v"^) 
2^ 



{n - 



2v{l - A) 



, (26) 
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where c is the normahzation factor and p(°)(n) is the 
Poisson distribution. By comparing Eq. H26I) with nu- 
merical simulations sX, v — \ (described below) we find 
that this formula accurately describes the actual form 
of P{n) in a much broader range of U /t than the solu- 
tion i(TT|i . 

D. Strong coupling limit 

Let us now turn to the strong coupling limit, vt/U 
1, at integer v, when atoms are well in the Mott in- 
sulator regime, and the zeroth order approximation to 
the wavefunction is a direct product of local Fock states, 
|*(t/;7)=o) — M^T-i — ^})- The effect of finite hopping t 
can be taken into account as a perturbation of the on-site 
interaction term in the Hamiltonian . To the first ap- 
proximation, this results in an admixture of the particle- 
hole pairs to \'^(t/u)=o), 

l*t/a> « (l + 4aj)\'^{t/u)=o), ^t/U « 1. 

<*,i> 

(27) 

With this wavefunction Eq. ^ leads to the following 
distribution 

P, = l-2P,_i, 
P.-i = Pu+i = 2dv{u + 1) t^/U^, 

Pn = 0, if \n-u\> 1. (28) 

E. Numerics 

Clearly, as U is changed from to cxd the number distri- 
bution must change qualitatively. Having an essentially 
long tail at 71 > 1/ in the weakly interacting limit, P„ 
becomes sharply peaked at large U with equal probabil- 
ities for — 1 and v + 1 particles on a site. To obtain 
Pn{U/t) in the crossover regime, t U, we perform MC 
simulations using the continuous-time Worm Algorithm 
scheme jl^. We set — 1 and take the limit Ns oo, 
[3 = 1/T — > oo, where T is the temperature. For the 

linear system size L = N^'^ = 24 and (3 ^ L/27rc, where 
c ~ 6t is the typical value of sound velocity in the super- 
fluid phase (higher temperatures can be used in the MI 
phase with gaps ~ U), the shape of the distribution is 
already well saturated, within a fraction of one percent 
in 3D and 2D, and within a few percent in ID, consistent 
with the fact that P„ is an essentially local characteristic. 

The simulation results for 3D, 2D and ID are shown 
in Fig.^ where Pq, Pi and P2 arc plotted as functions of 
U/t along with the predictions of Eq. and Eq. 
The main observation is that, although Pn{U/t) funda- 
mentally does not reveal any critical behavior, in 3D 
and 2D the SF-MI transition is marked by a substantial 
change in P„ curves — they rapidly plateau in the Mott 
regime for U/t > (U/t)c- As expected, in ID the curves 



are much more smooth and the saturation at high U is 
not pronounced. Remarkably, Eq. (|26|l . deviates notably 
(by a few percent) from the numerical data only at U/t 
as high as ~ 2 in ID, ~ 4 in 2D, and ~ 8 in 3D. 

We compare the theoretical predictions of section III Bl 
for the asymptotics in the v 1 limit with the results of 
Monte-Carlo simulations. From Fig. |21 where the sim- 
ulated probability distribution at [//t = 2 is plotted 
along with Gaussian curves with variances calculated by 
Eq. (|24|l at corresponding v, we conclude that the shape 
of the distribution is perfectly well (within the error bars) 
described by the Gaussian distribution already at v — b. 
In Fig.|31 the numerically simulated and the prediction 
of Eq. H24|l are plotted together as a function of U/t. 



III. COMPARISON WITH EXPERIMENT: 
EFFECTS OF FINITE TEMPERATURE 

In actual experiments, atoms are confined in the lat- 
tice by a trapping potential, typically of a parabolic form. 
This results in an inhomogeneous density profile in the SF 
regime and in the formation of Mott-plateaus 0,0,01 — 
spherical shells of integer filling — in deep lattices. Cor- 
respondingly, the distribution of the number fluctuations 
is also inhomogeneous, i.e. P^ = Pn{i), where i is the 
site index. The development of experimental techniques 
allowin g de tection of atoms with a single-site spatial res- 
olution [2g should open an exciting possibility to directly 
measure atom correlation functions and P„(«) in partic- 
ular. Current experiments |l. ll6llT7l | typically deal with 
integral characteristics of the number distribution, such 
as the fraction of the total number of atoms found on 
lattice sites with occupation n, 

/„-;^nP„(z). (29) 

i 

The fraction of pairs /2 in both SF and MI regimes can 
be accurately probed by the spin-changing collisions 0]- 
The measurement is set up in the following way [13. Af- 
ter the system is allowed to equilibrate at a fixed value of 
the lattice potential depth Vq, the configuration of atoms 
is frozen by a rapid increase of Vq. Then, coherent spin 
dynamics in the two-particle channel can be induced with 
a near-unitary efhciency, and the spin oscillation ampli- 
tude is measured to yield the fraction of pairs. With this 
technique, Gerbier et al. observed /2 as the system 
was driven from a SF regime across the transition deep 
into MI regime, corresponding to the values of the ini- 
tial lattice potential Vb ranging from AEr to 40£'r, where 
Er — h? /2m\^ is the single photon recoil energy, m is 
the mass of a ^^Rb atom, and A is the lattice laser wave- 
length. 

In our simulations, the system of ^^Rb atoms of 
Ref. ^3 is described by the Bose-Hubbard model iQ. 
The external potential is harmonic, = muj^x^ /2^ and 
all the parameters of the Hamiltonian J^l are fixed by 
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U/t 

FIG. 1: (Color online.) The probability P„ to detect n particles on a single site of a homogenous square lattice as a function 
of U/t at zero temperature and unitary filling. The results of MC simulation are shown for n — (squares), n = 1 (circles) 
and n — 2 (triangles). The uncertainties in P„ are smaller than the symbol size. The predictions of Eq. 12611 in the SF regime 
and that of Eq. 1281 in the strong coupling limit are plotted by thick and thin solid lines respectively. The critical points of 
the SF-MI transition, marked by the arrows on the graph, are (U/t)c = 29.34(2) in three dimensions ^23] (a), {U/t)c ~ 16.4(8) 
in two dimensions |3C| (b), and {U/t)c = 3.30(2) in one dimension [sj (c). 



0.1 







■ v=5 






• v=10 






A v=15 









u.u r ' ~^ ' — — r ' 1 — ' 1 — 

5 10 15 20 

Occupation number, n 

FIG. 2: (Color online.) Probability P„ of the on-site occupa- 
tion number n = rii + i/a.tU/t — 2 from the Monte-Carlo sim- 
ulation (the error bars are smaller then the symbol size) with 
the filling factor 5 (squares), 10 (circles), and 15 (triangles). 
The solid lines are the Gaussian curves with the variances 
calculated from Eq. 1241 with the parameters of simulations. 
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FIG. 3: (Color online.) The variance as a function of in- 
teraction strength U/t from the Monte-Carlo simulation (the 
error bars are smaller then the symbol size) with the filling 
factor 5 (squares), 10 (circles), and 15 (triangles). The solid 
lines represent the prediction of Eq. (1241 with the parameters 
of simulations. 



the experimental setup. We study the number fluctua- 
tions at two values of the lattice depth, Vq — 8Er and 
Vq = 13Er (corresponding to U/t ~ 7.4 and U/t « 35.6 
respectively), which serve as examples of typical SF and 
MI phases in the strongly correlated regime. The trap- 
ping frequency ujq is equal to 2t: x 37Hz and 2Tr x 46Hz for 
Vb = SEr and Vq = 13Er correspondingly p^l3^. We 



perform simulations in a sensible range of temperatures. 
Direct comparison between experimental and numerical 
data enables us to estimate the final temperature of the 
system. The results of the simulations are sown in Fig. 01 
where the curves for /i,/2, and /a as functions of the to- 
tal atom number in the trap N are parameterized by the 
temperature. 

Thermal effects vanish and a system can be considered 
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FIG. 4; (Color online.) Fractions of atoms /i, /2, /a occu- 
pying lattice sites with one, two and three particles on them 
respectively versus the total number of particles in the trap 
at the lattice depths Vb — 8Er, (a), (b), (c), and Vb = 13-Er, 
(d), (e), (f). The numerical results at different temperatures 
are plotted with the experimental data for /2 of Ref. ifi^j. 

to be in its ground state, if the temperature is substan- 
tially smaller than the energy of the low- lying excitations. 
In the SF regime, these modes have typical frequencies 
of order -^/m/m, uo [m* is the effective mass in the lat- 
tice), which for the parameters of Fig. ^a)-(c), gives a 
rough estimate T <C 3i. Clearly, the curves for T = t/i 
and T = t in Fig. Ela)-(c) are in the regime of effective 
zero temperature. Already at T = 2t the thermal effects 
become significant resulting in a larger atom cloud and, 
consequently, reduced average density. At T < 5t the 
size of the cloud for the maximum number of atoms in 
the trap is ~ 100 lattice sites in each dimension. To ob- 
tain /„ at r = lot, we have to resort to finite-size scaling, 
being limited by computer memory at the linear system 
size of 150. 

In the MI phase, the zero temperature regime is 
reached for temperatures smaller than the excitation gap, 
which is of order U/2 for the parameters of Fig.01d)-(f). 
This leads to the condition T ^ 20t, consistent with the 
saturation of /„ below T — 5t (see Fig. ^d)-(f)). For 
N < 5 X 10^, in the zero temperature limit, the curves 
/„ are flat, corresponding to the filling of the ly ~ 1 Mott 



shell, and the number distribution is essentially squeezed. 
Simple estimates QJj show that the decrease of /i ac- 
companied by the increase in /2 at iV '--^ 5 x 10^, and 
the saturation of /2 with a peak in /a at iV ~ 2 x lO'^ 
are consistent with the formation of Mott plateaus with 
ly = 2 and = 3, respectively. As seen from Fig.01[d),(e), 
final temperature effects degrade the degree of number 
squeezing in the MI by favoring particle-hole excitations. 

The comparison of the calculated fraction of atom pairs 
/2 with the measurements of Ref. ^3 (open circles in 
Fig. EIb),(e)) gives the typical experimental tempera- 
tures of order of 5t « 1.5 x 10~^ Er in the SF regime 
and lot « IQ-'^Er in the MI regime. Temperatures on 
the order of a few t have been also observed in a (one- 
dimensional) Tonks- Girardeau gas where the effec- 
tive fermionization due to strong interactions allows to 
deduce the experimental temperature from the momen- 
tum distribution. Note that the system acquires a finite 
temperature as a result of its loading into the optical lat- 
tice, and therefore the final temperature is supposed to 
depend on the number of atoms in the trap. However, 
the fact that the experimental data lie above the T = 
curve in Fig. ^b) at high N is unlikely to be explained 
by the the effects of heating. When a large /a fraction 
is present, a change in the spin resonance condition can 
result in a considerable contribution of spin collisions on 
triply o ccup ied sites to the observed spin oscillation am- 
plitude [13, which could explain the anomalously high 
apparent /2 [3^ . Such drifts of the resonance parame- 
ters are carefully checked for, but can not be ruled out 
completely [s^ . 



IV. CONCLUSIONS 

We studied the ground-state on-site number statistics 
of interacting lattice bosons. We considered the limits 
of weak interatomic interactions, the limit of large filling 
in the SF regime, and the limit of strong interatomic in- 
teractions. At lyU/t ^ 1, the correction to the Poisson 
distribution is described by Eq. (|ll|l . with an essentially 
dimension-dependent scaling — Eq. H15() in ID, Eq. H18() 
in 2D, and Eq. I|lt)|) in 3D. In the case of large occu- 
pation numbers, 3> 1, we show that, in the region 
of interactions U/t <C i^, the on-site occupation num- 
ber distribution is Gaussian and its variance, given by 
Eq. (|24|l . gradually decreases with the asymptotic form 
oc yJvTfU at l/i^ <C U/t <^ u hi all dimensions. An 
excellent agreement with numeric simulations is found 
already at ly = 5. At vt/U <C 1 and integer filling v, the 
distribution, Eq. (|28|l . is dominated by particle- hole fluc- 
tuations on top of an ideal MI. At = 1, we performed 
Monte Carlo simulations of the on-site occupation num- 
ber distribution in ID, 2D, and 3D in a wide range of U/t 
covering the crossover region between the limiting cases. 

We simulated a parabolically confined system in the 
realistic case of Ref. ^3 with the final temperature of 
the system being the only free parameter. By direct 
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comparison between experimental data and numerical re- 
sults at different temperatures, we were able to estimate 
the experimental temperature, which we found to be of 
the order of a few t near the transition. The error bars 
are small enough to determine T with accuracy of order 
t. Our results suggest that measurements of the on-site 
atom number fluctuations can serve as a reliable method 
of thermometry in both superfluid and Mott-insulator 
regimes. We believe that with more elaborate techniques, 



such as that giving access to the spatial number distri- 
bution Pn{i), the temperature resolution can be further 
improved. 

We are grateful to Fabrice Gerbier for providing us 
with the experimental parameters and for fruitful discus- 
sions of the results of this work. This research was sup- 
ported by the National Science Foundation under Grant 
No. PHY-0426881. 
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